function gipl = fgiplread(file,silent)
%
% GIPL reader for simple data
% MJ White, CMIC UCL <mark.white@ucl.ac.uk>
%
% Usage:
%   gipl = fgiplread(file)
%   gipl = fgiplread(file,silent)
%
% If silent is unspecified (and is zero), a couple of lines will be
% written to confirm the endianness and dimensions of the GIPL.
%
% Returns a GIPL structure with components:
%   gipl.dim
%   gipl.type
%   gipl.voxdim
%   gipl.desc    (returned as a string)
%   gipl.matrix  (returned as a 4*4 matrix)
%   gipl.imin
%   gipl.imax
%   gipl.origin
%   gipl.extended
%   gipl.magic
%
% For full description of possible GIPL header field meanings see
%   http://isen.cs.ucl.ac.uk/trac/wiki/FormatGipl
%
% All header fields (known to be in current use) are read and
% returned, but no particular validation is done because the actual
% range of use cases of the GIPL header fields is so wide.  Beware,
% for example, that imin and imax may not have been filled in when the
% file was written; and that matrix has several possible meanings.
%

  fmt = 'b';

  if (nargin < 2)
    silent = 0;
  end

  function n = mydisp(s)
    if silent == 0
      disp(s);
    end
    n = 0;
  end

  function v = check_magic(fh) 
    fseek(fh,252,'bof');
    magic = fread(fh,1,'uint32',0,'b');
    if magic == 719555000 | magic == 4026526128
      mydisp('Big-endian GIPL');
      v = 'b';
    elseif magic == 3096044330 | magic == 2968125423
      mydisp('Little-endian GIPL');
      v = 'l';
    else
      error 'Bad magic number, probably not a GIPL.'
    end
  end

  fh = fopen(file);
  gipl = struct();

  fmt = check_magic(fh);

  % read header block
  fseek(fh,0,'bof');
  gipl.dim      = fread(fh,4,'uint16',0,fmt);
  gipl.type     = fread(fh,1,'uint16',0,fmt);
  gipl.voxdim   = fread(fh,4,'float32',0,fmt);
  gipl.desc     = strcat(fread(fh,80,'*char',0,fmt));
  mtemp         = reshape( fread(fh,12,'float32',0,fmt), [4 3] )' ;
  gipl.matrix   = [ mtemp ; 0 0 0 1 ];
  fseek(fh,4,'cof');    % identifier
  fseek(fh,28,'cof');   % spare bytes
  fseek(fh,1,'cof');    % flag
  fseek(fh,1,'cof');    % orientation
  gipl.imin     = fread(fh,1,'float64',0,fmt);
  gipl.imax     = fread(fh,1,'float64',0,fmt);
  gipl.origin   = fread(fh,4,'float64',0,fmt);
  fseek(fh,4*2,'cof');  % unknown
  gipl.extended = fread(fh,1,'uint32',0,fmt);
  fseek(fh,4*1,'cof');  % unknown
  gipl.magic    = fread(fh,1,'uint32',0,fmt);

  % read data
  total = prod(gipl.dim);
  mydisp(sprintf(' Dimensions: %d %d %d %d', ...
    gipl.dim(1), gipl.dim(2), gipl.dim(3), gipl.dim(4)));
  cpx = 0;
  switch gipl.type
  case 1
    type = 'ubit1';
  case 7
    type = 'int8';
  case 8
    type = 'uint8';
  case 15
    type = 'int16';
  case 16
    type = 'uint16';
  case 31
    type = 'int32';
  case 32
    type = 'uint32';
  case 64
    type = 'float32';
  case 65
    type = 'float64';
  case 144
    type = 'int16';
    cpx = 1;
  case 160
    type = 'int32';
    cpx = 1;
  case 192
    type = 'float32';
    cpx = 1;
  case 193
    type = 'float64';
    cpx = 1;
  otherwise
    error('Unknown GIPL data type.');
  end
  if cpx == 0
    data = fread(fh,total,type,0,fmt);
    data = reshape(data,gipl.dim');
  else
    datar = fread(fh,total*2,type,0,fmt);
    datar = reshape(datar,[2 total]);
    data = complex(datar(1,:),datar(2,:));
    data = reshape(data,gipl.dim');
  end
  gipl.data = data;

  fclose(fh);
end
